Our MD simulations were performed using the graphics processing unit (GPU)-based MD package MDPSCU.^[31] For the inter-atomic potentials, we implemented the embedded-atom-method (EAM) potentials developed by Bonny et al.^[32] for the W–H–He system on the basis of the W–W potential of Marinica et al.^[33] Only the potential denoted as EAM1 in Ref. [32] was used. All the simulations were performed in the temperature range of 300–2000 K. For W, the lattice constant a₀ = 3.14 Å was adopted, which is consistent with that used by Marinica et al.^[33]

In calculating the H diffusion coefficient, 1000 independent simulation boxes were generated for body-centered-cubic (bcc) W. Each box has a size of 10a₀ × 10a₀ × 10a₀ and contains 2000 W atoms. Periodic boundary conditions were applied in the x-, y-, and z-directions. In the simulation procedure, we initially added one H atom to every simulation box, and then thermalized them to a given temperature by the Monte Carlo sampling of the Maxwell distribution for the velocities of the atoms. After achieving thermal equilibrium, the boxes were then allowed to freely relax for 40 ps with the integral time step of 0.25 fs. The time step is small enough to guarantee the energy conservation. Since the energies of the boxes follow the Boltzmann distribution at the given temperature, an NVT ensemble is simulated. The displacement R(t) of the H atom was recorded. The diffusion coefficient D at a given temperature was then calculated by linear fitting the equation 1 where t is the time and ⟨R²(t)⟩ denotes the ensemble average.

To study the diffusion behavior of H atoms near surfaces, we constructed W{100} and W{111} surfaces with bcc W box sizes of 10a₀ × 10a₀ × 30a₀ and 8.5a₀ × 9.8a₀ × 27.7a₀, respectively. Generally, the simulation procedure consisted of two phases. In the first phase, one H atom was randomly introduced into the constructed bcc W box, and periodic boundary conditions were used in all directions. For each of the box sizes, we generated 500 independent replicas that were subsequently thermalized to a given temperature. In the second phase, the periodic boundary condition in the z-direction was discarded, which is defined as the (100) and (111) crystal orientations. The length of the box side in the z-direction was sufficient to eliminate possible interactions between the two surfaces. In this phase, all the boxes were evolved for 1000 ps, during which the states of the boxes were recorded every 20 ps.

Using the methods described above, the diffusion coefficients D of H, D, and T in W were calculated in the temperature range of 300–2000 K and fitted to the Arrhenius relationship 2 as shown in Fig. 1. Here, D₀ is the pre-exponential factor (or pre-factor), E_a is the activation energy, k_B is the Boltzmann constant, and T is the temperature (K). At first sight, the dependence of D on temperature follows the Arrhenius relationship well, with the corresponding D₀ and E_a values for H, D, and T presented in Table 1. However, detailed analysis shows that E_a varies with the temperature range chosen when fitting the MD data to Eq. (2). For example, over the whole temperature range, fitting the MD data for an H atom to Eq. (2) results in E_a = 0.218±0.002 eV (Table 1), while the E_a values obtained in the discrete temperature ranges of 300–1000 K, 1000–2000 K, and 1600–2000 K are 0.221 ± 0.002 eV, 0.201 ± 0.004 eV, and 0.183±0.003 eV, respectively. In the low-temperature range (300–1000 K), the E_a value (0.221 eV) is in good agreement with the diffusion barrier of 0.221 eV that we obtained using the nudged elastic band (NEB) method.^[34,35] This can be understood because the NEB is a static method which calculates the lowest diffusion barrier on the potential surface at 0 K. However, with increasing temperature, E_a decreases, and could even be smaller than the diffusion barrier of the NEB. This variation of E_a is an indication of non-Arrhenius diffusion. A similar phenomenon has been observed for He in W by Wang et al.,^[25] Shu et al.,^[36] Perez et al.,^[27] and Wen et al.^[37] The E_a value of He atoms in W is lower in the temperature range 650–1000 K than that at 300–650 K.^[25,27,37] At temperatures higher than 1000 K, the diffusion of He atoms in W tends to be the Einstein–Smoluchowski type,^[37] which would result in an enlarged E_a if arbitrarily fitted to the Arrhenius relationship. The mechanism of multiple diffusion pathways, proposed by Wang et al.^[25] and Shu et al.,^[36] cannot explain the lowered E_a in the intermediate temperature range (650–1000 K). The reason is as follows. According to HTST for multiple diffusion pathways, the diffusion coefficient D can be expressed as 3 where D_i and E_i represent the pre-factor and activation energy of diffusion pathway i. E₀ is the lowest activation energy. The logarithm of diffusivity lnD can be expressed as 4 Because E₀ < E_i, lnD tends to be and at the limit T → 0. It is easy to prove that is always less than −E₀/k_B and monotonically decreases with increasing temperature T. In other words, the activation energy obtained by fitting the Arrhenius relationship would increase with increasing temperature. Wen and coworkers suggested that the deviation from HTST prediction could be the result of many-body effects.^[37] Alternatively, we conjecture that the existence of correlated successive jumps, which have been observed in the diffusion of metallic adatoms on surfaces,^[38] is a more likely explanation for the non-Arrhenius diffusion of He and H isotopes in W. In comparison with the non-Arrhenius diffusion of He in W, the non-Arrhenius behavior for H in W is insignificant, because the lowest diffusion barrier for H in W is higher than the lowest diffusion barrier for He in W. Non-Arrhenius diffusion for H in W tends to occur at a higher temperature than for He in W, and the transition to Einstein–Smoluchowski diffusion is not observed up to the highest temperature (2000 K) considered in the present paper. The non-Arrhenius diffusion of H in W has been more obviously observed in the MD simulations of Liu et al.,^[39] who used a different potential for the H–W system than ours. However, in contrast to our findings, they observed E_a increasing as temperature increased, and an E_a in the low temperature region that was less than the value they obtained by NEB using their potential.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) Diffusion coefficients of H, D, and T in W, calculated by MD simulation (DH, DD, and DT), by fitting the Arrhenius relationship to the MD data (lines). is the diffusion coefficient of T, derived from DH based on single-mode harmonic transition state theory.

Table 1.

Table 1.

Table 1. Fitted pre-exponential factors and activation energies for H, D, and T and a comparison with literature data. . H isotopes D0/m2·s−1 Ea/eV NEB: Ea/eV H (8.03 ± 0.012) × 10−8 0.218 ± 0.002 0.211 Present work D (7.551 ± 0.025) × 10−8 0.225 ± 0.001 0.211 T (6.902 ± 0.017) × 10−8 0.227 ± 0.001 0.211 Frauenfelder[12] H 4.1 × 10−7 0.391 (1100–2500 K) Frauenfelder[22] H 1.58 × 10−7 0.25 (1500–2500 K) Heinola et al.[22] H 5.2 × 10−8 0.21 Johnson et al.[21] H 8.93 × 10−7 0.38

Table 1. Fitted pre-exponential factors and activation energies for H, D, and T and a comparison with literature data. .

For comparison with experimental results, the activation energies obtained by separately fitting the Arrhenius relationship to Frauenfelder’s experimental data in the temperature ranges of 1100–2500 K and 1500–2500 K are 0.39 eV and 0.25 eV, respectively.^[22] The observation that the E_a at higher temperatures is smaller than that at lower temperatures is similar to our results shown above. However, the effect of non-Arrhenius diffusion on the experimental observations remains ambiguous because of the large discrepancies between E_a values obtained through experiments and simulations.

According to the HTST of Vineyard,^[40] E_a is independent of isotopic mass, and the pre-factor in Eq. (2) for isotopes can be written as 5 where α denotes the isotopes (H, T, D), A is a constant, and is the effective mass. Vineyard showed that is dependent on the masses of the jumping atom (here, H, D, or T) and the host atom (W), as well as the direction of the diffusion path at the saddle point on the potential surface. would be well approximated by the mass of the jumping atom if its movement is the dominant contribution to the configuration variation during jumping. By this approximation, the diffusion coefficients for D and T can be expressed by 6 where D_H is the diffusion coefficient of the H atom. In Fig. 1, the diffusion coefficient is compared to the diffusion coefficient D_T obtained by MD simulation. agrees well with D_T at low temperature but deviates from D_T with increasing temperature. The difference arises from two sources, as discussed below.

First, from Table 1, the pre-factors calculated via the MD simulations give the ratio D_0,T/D_0,H = 0.93. This value is larger than , the ratio predicted by Eq. (6). From Eq. (5), when the value of D_0,T/D_{0, H} approaches unity, the effective masses and are similar. According to Vineyard’s analysis,^[40] the underlying physical mechanism is that many-body effects play roles: the W atoms neighboring the jumping atom (H or T) are displaced non-symmetrically as the jumping atom passes through the saddle point of the potential surface, and thus, the W atoms impart a large contribution to the effective mass. The situation is similar for the diffusion of D, with D_{0, D}/D_{0, H} = 0.97 vs. . Second, also from Table 1, differences exist among the E_a values obtained by MD simulations for H, D, and T. Although the differences are small, they contradict the predictions of the HTST, in which E_a is constant for isotopes. Since MD is based on classical mechanics, the differences among the E_a values of the isotopes do not arise from the contribution of the zero-point energy (ZPE), which is an isotope-dependent quantum effect that was added to E_a in previous studies.^[21,22] Moreover, it has been shown that including the ZPE in E_a leads to a more temperature-dependent diffusion-coefficient pre-factor than when the ZPE is not included.^[23] Clearly, the dependence of E_a on isotope identity observed in the present paper derives from the anharmonic properties of the potential surface.

MD simulations were performed for temperatures from 600–1600 K with the simulation boxes having two free surfaces. Two surface orientations, {100} and {111}, were considered.

For the case of the W{100} surface, figure 2(a) displays snapshots taken at times t = 0 and 1000 ps for temperatures T = 1000 K and 1400 K. The snapshots were generated by merging snapshots of 500 simulation boxes. From the figure, the interstitial H atoms are initially uniformly distributed in the simulation boxes. As time elapses, some of the H atoms escape from the substrate, and dilute H atoms are found in the subsurface region. Figure 2(b) shows the corresponding depth distribution of H atoms. These results suggest that H atoms near the surface can quickly escape from the substrate, analogous to the behavior of He atoms near the W{001} surface.^[30] The time scale for the escape of subsurface H atoms should be shorter than the recording interval (20 ps).

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) (a) Snapshots of H-atom distribution in W bulk for W(100); t = 0 and 1000 ps; and T = 1000 K and 1400 K. Each picture represents the merging of 500 simulation boxes. The red and green spheres represent W and H atoms, respectively. (b) The corresponding depth distribution of H atoms.

Figure 3 displays the snapshots and depth distribution of H atoms for the case of the W{111} surface orientation, which is in the closest-packed direction of W atoms. In the case of He atoms in W, temperature-dependent trap mutations may occur when He atoms move to the W{111} surface, and these trap mutations induce the subsurface accumulation of He atoms.^[28,30] We repeated the MD simulation of He diffusion near the W{111} surface using the W–H–He potential of Bonny et al.,^[32] which is used in the present paper for the W–H interaction but is different from the W–He potentials used in Ref. [28] and [30]. The simulation results (not shown in the present paper) further verify the accumulation of He atoms in the subsurface region of the W{111} surface. However, in contrast to the He case, no subsurface accumulation of H atoms is observed in Fig. 3. Similarly to the case of H diffusion near the W{100} surface described above, H atoms can quickly escape from the substrate through the W{111} surface. This observation is applicable to all temperatures considered in the present paper.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) (a) Snapshots of H-atom distribution in W bulk for W(111); t = 0 and 1000 ps; and T = 1000 K and 1400 K. Each picture represents the merging of 500 simulation boxes. The red and green spheres represent W and H atoms, respectively. (b) The corresponding depth distribution of H atoms.

Because the time scale for the escape of subsurface H atoms from the substrate is much smaller than the time scale considered in continuous-diffusion-theory (CDT)-based simulations (e.g., the size of the time grid in rate theory), we can establish a model for continuous diffusion for W surfaces. In this model, the W surface is treated as an absorber. The H atoms can be considered to migrate through the bulk tungsten until arriving at the absorber, where they are absorbed instantly with probability R. According to Wang et al.,^[30] R can be determined by fitting the MD data and the following formula 7 with where N_a(t) is the accumulated number of H atoms escaping from the substrate in time t, and N_B is the number of simulation boxes, which equals the total H atoms (here, N_B = 500). D is the diffusion coefficient of H in the W bulk calculated above. The n₀(z) term represents the initial depth distribution of the H atoms and is given by 1/L_z (L_z is the simulation box size in the z-direction) since the H atoms are initially uniformly distributed in the simulation boxes. h is the distance from the center of the simulation box to the absorbing layer, and is defined as h = L_z/2 here.

The absorbing probabilities R determined for the W{100} and W{111} surfaces are summarized in Table 2 for different substrate temperatures. Figure 4 shows N_a(t)/N_B obtained by MD simulations as well as by Eq. (7). Generally, the MD results are well described by Eq. (7), suggesting validation of the absorber model. By CDT, the value of R should be 1 for a perfect absorber. As we can see from Table 2, the R value for W{100} is less than 1 while that for W{111} is greater than 1. The deviation of the R value from 1 reveals that the diffusion barriers of H atoms near the surfaces vary relative to the in-bulk diffusion barriers. Using the NEB method, we calculated the potential evolution for an H atom moving from the near-surface region to the top of the surface. The results are shown in Fig. 5. It should be noted that the horizontal axis represents the configuration images of the NEB rather than the depth of the H atoms in W. Thus, for the W{111} surface, the interstitial state of an H atom close to the surface has a potential energy lower than that of the in-bulk interstitial state by 0.07 eV, while for the W{100} surface, the interstitial state of a near-surface H atom has almost the same potential energy as the in-bulk interstitial state. This observation explains why the R value is larger than 1 for the W{111} surface and smaller than 1 for the W{100} surface. Meanwhile, the energy difference from the bulk to the surface is approximately 1.50 eV, which is in good agreement with the absorption energy (from the surface to the bulk) calculated by Hodille et al.^[41]

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) The ratio of escaping H atoms (Na) to the total H atoms (NB) as a function of time obtained from MD simulations for (a) W{100} and (b) W{111} surfaces. The solid lines are the curves generated from the fitting of Eq. (10) to the MD data.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Evolution of the potential barrier for H-atom movement in the near-surface region until escaping from the surface: (a) for W{100} and (b) for W{111} surfaces. The potential energy of the in-bulk tetrahedral interstitial state is taken as zero.

Table 2.

Table 2.

Table 2. R values obtained by fitting Eq. (10) to simulation data. The diffusion coefficient used in Eq. (10) is calculated in the present work. . Temperature/K Diffusion coefficient/m2·s−1 R for W{100} R for W{111} 600 1.244 × 10−9 0.809 0.937 800 3.607 × 10−9 0.944 1.232 1000 6.593 × 10−9 0.925 1.140 1200 9.796 × 10−9 0.977 1.086 1400 1.347 × 10−8 0.938 1.097 1600 1.616 × 10−8 0.931 1.112

Table 2. R values obtained by fitting Eq. (10) to simulation data. The diffusion coefficient used in Eq. (10) is calculated in the present work. .

We note that the W–H potential used in this paper is repulsive if H atoms are on the top of the W surfaces. No H adatoms, and thus, no H₂ molecules, which are thought to be detected in TDS experiments, are observed on the top of the W surfaces in our MD simulations. The study of the diffusion, recombination, and desorption of H adatoms on W surfaces likely needs a H–W potential that can correctly describe the interaction between H adatoms and W surface atoms. However, the model established here can be applied to modelling the production rate of H adatoms on W surfaces produced by diffusing in-bulk H atoms. The behavior of H adatoms on W surfaces can be considered as an independent subsequent process and beyond the scope of the present paper.